#include "TestLinear.h"
#include "RungeKutta.h"
#include <cmath>
#include <iostream>
#include <cstdio>

using namespace std;

const char* TestLinear::DOUBLE_SINGLETON = "%-20.16f\n";
const char* TestLinear::DOUBLE_TRIPLET = "%-20.16f%-20.16f%-20.16f\n";

double TestLinear::f1(const double* x) {
	return 2 * x[1] + x[2] + x[3];
}

double TestLinear::f2(const double* x) {
	return -2 * x[1] - x[3];
}

double TestLinear::f3(const double* x) {
	return 2 * x[1] + x[2] + 2 * x[3];
}

void TestLinear::run() {
	const int dim = 3;
	const double a = 0;
	const double b = 1;
	const double h = 1E-4;

	vector<int> stepsCounts;
	stepsCounts.push_back(10000);
	Mesh mesh{ a, b, &stepsCounts };

	MultivarFunction rightHands[] = {
		[this](double* x) {
			return f1(x);
		},

		[this](double* x) {
			return f2(x);
		},

		[this](double* x) {
			return f3(x);
		},
	};

	double initialValues[] = { 0, -1, 2 };

	RungeKutta rungeKutta;
	rungeKutta.setDim(dim);
	rungeKutta.setMesh(&mesh);
	rungeKutta.setRightHands(rightHands);
	rungeKutta.setInitialValues(initialValues);
	
	vector<NodeSolution*>* y = rungeKutta.solve();

	auto u1 = [](double x) {
		return -exp(x) + exp(2 * x);
	};

	auto u2 = [](double x) {
		return (1 + x) * exp(x) - 2 * exp(2 * x);
	};

	auto u3 = [](double x) {
		return -x * exp(x) + 2 * exp(2 * x);
	};

	for (const auto& nodeSolution : *y) {
		double x = nodeSolution->getX();
		printf(DOUBLE_SINGLETON, x);

		double* yi = nodeSolution->getY();
		printf(DOUBLE_TRIPLET, yi[0], yi[1], yi[2]);
		
		printf(DOUBLE_TRIPLET, u1(x), u2(x), u3(x));
		
		cout << '\n';
	}

	Utils::freeSolution(y);
}